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Vulnerability Analysis of Power Systems 

Taedong Kim' , Stephen J. WrighP, Daniel Bienstock*, and Sean Harnett !i 


Abstract —Potential vulnerabilities in a power grid can be 
exposed by identifying those transmission lines on which attacks 
(in the form of interference with their transmission capabilities) 
causes maximum disruption to the grid. In this study, we 
model the grid by (nonlinear) AC power flow equations, and 
assume that attacks take the form of increased impedance along 
transmission lines. We quantify disruption in several different 
ways, including (a) overall deviation of the voltages at the 
buses from 1.0 per unit (p.u.), and (b) the minimal amount of 
load that must be shed in order to restore the grid to stable 
operation. We describe optimization formulations of the problem 
of finding the most disruptive attack, which are either nonlinear 
programing problems or nonlinear bilevel optimization problems, 
and describe customized algorithms for solving these problems. 
Experimental results on the IEEE 118-Bus system and a Polish 
2383-Bus system are presented. 

Index Terms —AC power flow equations, vulnerability analysis, 
transmission line attack, bilevel optimization. 

I. Introduction 

Identifying the vulnerable components in a power grid is 
vital to the design and operation of a secure, stable system. 
One aspect of vulnerability analysis is to identify those trans¬ 
mission lines for which minor perturbations in their conductive 
properties leads to major disruptions to the grid, such as 
voltage drops, or the need for load shedding at demand nodes 
to restore feasible operation. 

Vulnerability assessment for power systems has been widely 
studied in recent times. Most works focus on minimizing 
the costs of load shedding and additional generation in the 
DC model (which is relatively easy to solve) or in lossless 
AC models (still relatively easy to solve and analyze). In 
ID, 0. 0, identification of critical components of a power 
system is formulated in a mixed-integer bilevel programming 
framework, and attacks on different types of system compo¬ 
nents (transmission lines, generators, and transformers) are 
considered. The lower-level problem in the bilevel formulation 
is replaced by its dual in |2j and is approximated using KKT 
conditions in 0. As an extension of 0, an approach based on 
Bender’s decomposition is proposed to solve larger instances 
of the transmission line attack problem in ||4). 

Vulnerability assessment using the lossless AC model is 
studied in 0, 0- In these papers, transmission-line attacks 
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are formulated as bilevel optimization problems, in which 
either unmet demands are maximized or attack costs (number 
of lines to attack) are minimized to meet a specified level of 
grid disruption. (These models are also discussed in 0, which 
describes the equivalence of the two models.) Then the lower- 
level problem is replaced by its KKT conditions, yielding the 
single-level optimization problem that is actually solved. In 
0, this mixed-integer problem is relaxed to a continuous 
problem (the binary variables are relaxed to real variables 
confined to the interval [0,1]), while 0 develops a graph¬ 
partitioning approach by identifying load-rich and generation- 
rich regions. 

The paper 0 describes a model that uses both load shed¬ 
ding and line switching as defensive operations to reduce the 
disruption of the system; the model is solved via Bender’s 
decomposition with a restart framework. Use of a genetic 
algorithm to solve the “N — k” problem (identifying the set of 
k lines in a grid of N lines whose removal causes maximum 
disruption) is discussed in 0. A minimum-cardinality ap¬ 
proach (solved using a cutting-plane method) and a continuous 
nonlinear attack model employing the DC power flow to 
represent power grids, where a fictitious adversary modifies 
reactances, are applied to the “N — k” problem in ma¬ 
in this paper, we propose two optimization models for 
vulnerability analysis. Both models are founded on the AC 
power flow equations, and both consider attacks in which 
the impedances of transmission lines are increased. In both 
formulations, the attacks respect a certain “budget;” their total 
amount of impedance adjustment cannot exceed a certain spec¬ 
ified level. The goal of the attacks is to maximize disruption, 
as measured by two different metrics. 

The first metric quantifies voltage disturbance at the buses, 
leading to a nonlinear programming formulation. The voltage 
disturbance usually appears as voltage drop, which often leads 
to an undesirable situation where voltages become low enough 
that the system cannot maintain stability. This situation, which 
is called voltage collapse or voltage instability, can happen 
either quickly or relatively slowly, and is characterized by a 
parallel process where reactive power demand correspondingly 
increases. This eventuality causes the active-power behavior 
of the system to approach the ’’nose” of the P — V curve. A 
more complete description is provided in ED pages 31 and 
35]. With our first metric, we estimate this possible voltage 
instability of a power grid assuming there is no response from 
a system operator to the attack. 

The second metric we consider here is a weighted sum of 
the amount of load shedding (at demand nodes) and generation 
reduction (at generation nodes) that is required to restore 
feasible operation of the grid following the attack. This power 
adjustment is considered as a defensive action of a system 
operator to keep voltages within a stable range to avoid 
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voltage collapse. This case is modeled as a bilevel optimization 
problem, in which the lower level finds the minimum load 
adjustment required to respond to the attack, and the upper- 
level problem is to find the most disruptive attack. 

In some existing literature, including some of the papers 
cited above, a bilevel optimization model is reformulated as a 
single-level optimization problem by replacing the lower level 
problem by its optimality conditions. This formulation strategy 
is unappealing, as the optimality conditions characterize only 
a stationary point, rather than a minimizer, so they may 
allow consideration of saddle points or local maximizers. In 
addition, if the bilevel formulation is designed to solve the 
attacker-defender framework that we consider in this paper, 
the reformulated single-level optimization model constructed 
by replacing the lower-level problem by primal-dual optimality 
conditions has the serious flaw that the model may exclude the 
most effective attack. Specifically, an attack (upper level deci¬ 
sion) that leads to an infeasible lower level problem obviously 
maximizes the disruption and thus is “optimal” for the attack 
problem (since it is not possible to make a operational decision 
at the lower level to defend against the attack). However, 
such an attack is excluded from consideration by the single- 
level reformulation since no (lower-level) primal-dual point 
satisfies the optimality condition constraints under the attack. 
Thus, the single-level formulation will ignore the most critical 
attack. Another drawback of single-level reformulations is 
that the optimality-condition constraints may violate constraint 
qualifications, causing possible complications in convergence 
behavior. 

The main contributions of this paper can be summarized as 
follows: 


1) In contrast to previous attack models, the grid is modeled 
with full AC power flow equations, which are the most 
accurate mathematical models of power flow. 

2) In our bilevel optimization formulation, we actually 
solve the lower-level problem rather than replacing it 
by its optimality conditions, as is done in earlier works, 
to avoid the formulation defects discussed above. 

3) We develop effective heuristics that make our formula¬ 
tions tractable even for power grids with thousands of 
buses. 


The remaining sections are organized as follows. We de¬ 
velop the optimization models in Section [11] and describe 


the challenges to be addressed in solving them. Section III 


describes heuristics and optimization techniques that address 
these challenges and that yield solutions of the problems. 
Experimental results on 118-Bus and 2383-Bus cases are pre¬ 
sented in Section[|V] and we discuss conclusions in Section[V] 


II. Problem Description 


In this section we discuss power systems background and 
notation, and describe our two formulations of the vulnera¬ 
bility analysis problem. We describe notation and background 
on power flow equations in Subsection |II-A| Our first vul¬ 
nerability model, based on a voltage disturbance objective, is 
discussed in Subsection |II-B| The second model, based on a 


power-adjustment criterion, is discussed in Subsection II-C 


A. Notations and Background 


We summarize here the power systems notation used in later 
sections, most of which is standard. 

• Set of buses: AT 

• Set of generators: Q c Af 

• Set of demand buses: V C Af 

• Index of the slack bus: s £ AT 

• Set of transmission lines: £ C Af x Af 

• Unit imaginary number: j = y/—l 

• Complex power at bus i £ Af: Pi + ]Q, (active power: P,: 

real power: Qf) 

• Complex voltage at bus i £ Af: Vie-’ 6i (voltage magnitude: 

Vj; phase angle: 9i ) 

• Difference of angles 9i and 9^, for (i,i r ) £ £: 9w := 

9i - 9i> 

• (*, i') entry of the admittance matrix for the (unperturbed) 

grid: Gw + jBw (conductance: Gw', susceptance: Bw). 
We assume that the set of generators Q and the set of demand 
buses V form a partition of Af. 

An attack on the grid is specified by means of a line 

\C.\ 

perturbation vector: 7 £ K+ , with 7 w denoting the relative 
increase in impedance on line (i, i') £ £. Specifically, an 
attack designated by the vector 7 causes conductances and 
susceptances to be modified as follows: 


Gw( 7) 


Bw(. 7) 



where /ij^j is the shunt (line charging) admittance of line 
( i,m ) £ £. More details on the bus admittance matrix can 
be found from lfl2l Chapter 9]. Note that when 7 w = 0 for 
all ( i , 1 ) £ £, the conductances and susceptances all attain 
their original (unperturbed) values. 

The AC power flow equations with perturbations 7 can be 
written as follows: 


F p {v,e- 7 )- 
F Q {V,9' 7 )_ 


= 0 , 


( 1 ) 


where the entries of F p and F () - (for all i £ Af) are defined 
as follows, for all i £ Af: 

F z P (V,9', 7 ) := 

Vi T. Vi' {Gw ( 7 ) cos + B w ( 7 ) sin 0 W ) - Pi, ( 2a ) 

i' 

if (V,0; 7 ) := 

Vi ^2 V 1 (Gw ( 7 ) sin 9w - Bw ( 7 ) cos 9w) ~ Qi ■ 

i' 


We assume throughout the paper that P t > 0 for generator 
buses i £ Q and Pi < 0 and Qi < 0 for demand buses 
i £ V. The power flow problem is to find the values of 












3 


the vectors V, 9 , P , and Q that satisfy equations (|2j, given 
the load demands I’d and Qd at load buses and the voltage 
magnitudes Vg and active power injection Pg at the generator 
buses. Conventionally, the reactive powers Qg are eliminated 
from the problem (since they can be obtained explicitly from 
( | 2 b| > for i £ Q, and appear in no other equations), yielding the 
following reduced formulation: 


7) 


-F^v^y 

f£(V,0; 7 ) 

-Fg(v,0-,y). 


= 0 . 


( 3 ) 


Here, V s , 9 S , Vg, I’g, Pd, and Qd are parameters associated 
with the network; 7 is the impedance modification vector 
described above; and Vd and 9 guv are the variables in the 
model. These equations usually can be solved using Newton’s 
method, when the system has a solution. For additional details 
of formulation of power flow problems, see fl2l Chapter 10]. 


II. Voltage Disturbance Model 

The AC power flow problem ([3]) often has multiple solutions 
ED. but only those solutions with Vi ~ 1.0 per unit (p.u.) 
for all i £ V are stable and operational in practice. In the 
vulnerability model described in this subsection, we use the 
sum-of-squares deviation Ty of the voltages from 1.0 p.u. as 
a measure of the disruption caused by an attack: 


■TV (7) “ 



( 4 ) 


.2 where V is obtained by 
solving F{V,9\y) = 0, 

when F(V, 9; 7 ) = 0 has 

no solution. 

Here, Ty is a function of 7, the vector of relative impedance 
increases. Note that only the voltage magnitudes of demand 
buses V are considered in Fy(y), since the voltage magni¬ 
tudes for generators and slack bus are given and fixed. We set 
•TV ( 7 ) = +°° when the attack results in an infeasible grid, 
since such attacks are the best possible. 

To limit the power of the purported attacker, we impose a 
constraint on the vector 7 , and define the voltage disturbance 
vulnerability problem as follows: 


Hv(k, 7 ) := max Fv{l) (5a) 

7 

subject to e T 7 < ny (5b) 

0 < 7 < ye, (5c) 

where e = ( 1 , 1 ,..., 1 ) T , the scalar 7 is an upper bound 
on relative impedance perturbation for each line, and k is 
the maximum number of lines that can be attacked at the 
maximum level. (Note that the actual number of lines attacked 
may be greater than n if non-maximal attacks are made on 
some lines.) 


Note that although the following model is 
alternative to ([5]), it is actually not valid: 

a plausible 

v ™ ax 

v-d,0x>uS,7 * — 

( 6 a) 

subject to F(V, 9- 7 ) = 0 

( 6 b) 

T 

e 7 < ny 

( 6 c) 

0 < 7 < ye. 

( 6 d) 


The reason is that when there is an attack 7 satisfying ( |5b > 
and ( |5c] ) that results in an infeasible grid, the formulation ([5 1 
will find it (with an objective function of +00) while the 
formulation © will not. In other words, the formulation © 
does not fully capture the adversarial nature of the attack. 
However, as a practical matter, these two formulations find the 
same solution in cases where every 7 satisfying ( [5b] ) and ( [5c] ) 
allows for a feasible solution of the AC power flow equations. 


C. Power-Adjustment Model 

Our second way to measure severity of an attack is to 
consider the minimum adjustments to power that must be made 
to restore the grid to feasible operation. Power adjustments 
take the form of shedding load at demand nodes and adjusting 
generation at generator nodes. (We use weights in the objective 
to discourage adjustment on nodes where it is undesirable, 
such as at generators whose output cannot be adjusted or 
at critical demand nodes whose load cannot be changed.) 
Calculation of this weighted sum of power adjustments in¬ 
volves solving a nonlinear programming problem that we call 
the feasibility restoration problem. This problem forms the 
lower-level problem in the bilevel optimization problem, as 
we outline at the end of this subsection. 

1) Feasibility Restoration: When the attack represented 
by 7 is too severe, the AC power flow equations 0 may 
not have a solution for which the voltages lie within an 
acceptable range. The feasibility restoration problem finds 
minimal adjustments to the power demands (at demand nodes 
V) and power generation (at generator nodes Q) for which 
feasibility is restored to the AC power flow equations. The 
formulation is as follows: 


FA 7 ) := min 

Vd ,6t>u q j 

G q g iPD) 

subject to 


+ a i ) + ^Z^i\ p i\Pi 

ieS iev 

p g {V, 9\y) — \Pg\O (erg — &g) = 
F£(V,9-,y)-\P v \Qp v = 0 
F%(V,9-,y)-\Qv\&Pv=0 
V < V D < V 


0 < ct+ < a+ 

0 < Ug < Cfg 

0 < PD < Ptd 


(7a) 

0 (7b) 
(7c) 
(7d) 
(7e) 
(7f) 
(7g) 
(7h) 


where aQb is element-wise multiplication of vectors a and b. 
Here, the variables a + ,a~, and p represent relative changes in 
demand loads and power generation, so that constraints ( |7b| , 
© and ( [7d| represent power flow equations |3]i in which the 
loads Pg, P-d, and Qd are modified. The parameters uj, repre¬ 
sent positive weights on the changes to loads and generation, 
indicating the desirability or undesirability of changes to that 
node. We note the following points. 

• The same variable pt is used in the active and reactive 
power balance equations © and m since active and 
reactive load shedding should occur in the same fraction. 

• Bound constraints ( |7f| >, ( f7g] i, and ( |7h| > on the load shed¬ 
ding variables limit the adjustments to a reasonable range 
(which may be zero for some buses). 
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• The weights w, could be set to large positive values to 
discourage changes on that node, and to smaller values 
when change is acceptable. The case in which no change 
at all is allowable on that node can be handled by setting 
the upper bound to zero in ( |7f| ), ( |7g| ), or ©■ Throughout 
the paper, we assume that cu, = 1 for all i, but note that 
other positive values of these weights can be used without 
any complication to the model. 

• Power generation at the generator nodes may be either 
increased or decreased in general, but the loads at demand 
nodes can only decrease. (Upper bounds af, a~, and ~p i 
should not exceed 1. This means that the type of a bus 
— generator or demand bus — cannot be changed.) 

• The bounds © guarantee that voltage levels are opera¬ 
tionally viable. 

The objective to be minimized in ([7]) is the weighted sum of 
power adjustments that are necessary to restore feasibility to 
the power flow equations. We define Fl(i) = +oo when it 
is not possible to restore feasibility by adjusting loads and 
generations (which usually happens because the constraints 
regarding acceptable voltage levels © cannot be satisfied 
even when load shedding is allowed). 

The feasibility restoration problem 0 is a nonconvex 
smooth constrained optimization problem in general, so we 
can expect to find only a local solution when using standard 
algorithms for such problems. The problem generalizes 0 in 
that if a solution of the latter problem exists, it will yield a 
global solution of (|7]i with an objective of zero when we set 
af = of = 0 for i 6 5 and pi = 0 for * £ V, provided the 
voltage constraints © are satisfied. Moreover, by the well- 
known sparsity property induced by l\ objectives, we expect 
few of the components of Og, Og, and px> to be nonzero 
at a typical solution of 0. The problem 0 may also have 
operational relevance, guiding the grid operator toward a set 
of decisions that can restore stable operation of the grid with 
minimal disruption. 

For convenience of later discussion, we state 0 in the 
following more compact form: 


F L {l) '■= min p T y ( 8 a) 

x,y 

subject to Fl(x, y; 7 ) = 0 ( 8 b) 

x<x<x ( 8 c) 

0 < y < y, ( 8 d) 


where F^(x, y, 7 ) = 0 represents the equality constraints ©- 
( |7dl l, x includes the circuit variables V and 9, and y includes 
the power-adjustment variables o + , cr ~, and p. 

2) Bilevel Formulation: The bilevel optimization formu¬ 
lation seeks the attack 7 for which the power-adjustment 
objective 7) is maximized subject to the same attack budget 
constraints as in (|5j, that is. 


%l(k, 7 ) := max 

T 

Mi) 

(9a) 

subject to 

e T 7 < «7 

(9b) 


0 < 7 < fe. 

(9c) 


By substituting from (| 8 j, we obtain a max-min problem: 


= max min 

7 x,y 

T 

P y 

( 10 a) 

subject to 

0 

II 

£ 

( 10 b) 


X < X < X 

( 10 c) 


0 <y<y 

( 10 d) 


e T 7 < nf 

( 10 e) 


0 < 7 < ie. 

( 10 f) 


Bilevel optimization problems are, in general, difficult to 
solve. For problems of the form ( fT0| i, it is possible for the 
upper-level objective Fl to change discontinuously at some 
values of 7 , even when the constraint function F k is smooth 
and nonlinear. 

For the power-adjustment formulation, there is an additional 
complication: For most feasible values of 7 , the objective 
function is zero. This is because power grids are often robust 
to small perturbations, so when even when many impedances 
change, it is often possible to continue meeting all demands 
while respecting operational limits on the voltage values. This 
feature makes it difficult to search for the optimal 7 , since 
it is difficult even to find a starting value of 7 that causes 
nonzero disruption. We have developed specialized heuristics 
to address this issue; these are described in Section |1II-C| 

III. Algorithm Description 

We discuss a first-order method for the following formula¬ 
tion, which generalizes 0 and 0 : 


H(k, 7 ) := max 

7 

Hi) 

( 11 a) 

subject to 

e T 7 < k 7 

(lib) 


0 < 7 < 7 e. 

( 11 c) 


Although the objective F is not convex or smooth, we solve 
it with the classical Frank-Wolfe method (also known as the 
conditional gradient method), which we describe in the next 
subsection. 

A. Frank-Wolfe Algorithm 

The Frank-Wolfe algorithm lfl4l solves a sequence of sub¬ 
problems in which a first-order approximation to the objective 
around the current iterate is minimized over the given feasible 
set. If the objective F in CD were smooth, we would solve 
the following problem at the fcth iterate 'y k : 

w k := argmax (g k ) T (w- i k ) 

“ (12) 
subject to e T w < tvy, 0 < w < 7, 

where g k is a gradient F( 7 ) at 'y k . The new iterate is obtained 
by setting 

7 fc+ 1 = l k + a k {w k - 7 fc ), 

for some a k £ (0,1]. (Frank and Wolfe lfl4l give a specific 
formula for a k that guarantees a sublinear convergence rate for 
smooth convex F. An exact line search would yield a similar 
rate.) Because of the special nature of our constraint set, the 
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problem ( fl2| ) is a linear program with a closed-form solution, 
whose components w k are defined as follows: 

k J 7 if g k is one of k largest entries in g k , 

1 0 otherwise. 

We determine the step size ak by a standard backtracking 
procedure. Given a constant (b £ (0,1), and starting from a = 
1 , we decrease the step size by a f- (bo until the following 
sufficient decrease condition is satisfied for a small Ci £ (0,1). 

T(y k + a{w k - 7 fe )) > T(y k ) + ci ag kT (w k - 7 fe ). (14) 

We define o k to be the value of a accepted by this criterion. 
The algorithmic framework is shown in Algorithm [T] We 
terminate when the step ( w k — 7 fc ) becomes small, or when 
the step size a becomes less than a predefined > 0. 

Convergence behavior of the Frank-Wolfe procedure with 
backtracking line search for the smooth nonconvex case has 
been analyzed by Dunn lfl5l Theorem 4.1], where it is shown 
that accumulation points are stationary. (This result does not 
apply directly to our cases, because of potential nonsmooth¬ 
ness of the objectives.) 


Algorithm 1 Vulnerability Analysis 

Require: 

7 : Upper bound for impedance increases 7 ,, i £ C\ 
k\ Number of lines to attack; 

7 0 : Feasible initial value of 7 ; 

Ensure: 

7 *: Impedance vector that optimizes the attack; 


1 

2 

3 

4 

5 

6 

7 

8 
9 


k <- 0; 

while k <MaxIter do 

Find the gradient g k of the objective T at ; 

Find linearized optimum w k from 0 

Use the backtracking to find step size a & £ [0,1]; 

7 fc+1 ^7 fe +afe(^ fe -7 fe ); 
k 4 — k 1; 

Stop if termination conditions are satisfied, and set 7 * •<— 7 fc ; 

end while 


B. Gradient Calculation 

Algorithm [l] requires calculation of a gradient g k of the 
objective function T at the current iterate 7 fc . We have noted 
already that the power-adjustment objective T = Tl may be 
nonsmooth, due to changes in the active set of the subproblem 
so the gradient may not be well defined. We note however 
that bF[, can reasonably be assumed to be smooth almost 
everywhere; changes to the active set can be expected to 
happen only on a set of measure zero in the feasible space 
for 7 . Our algorithm does not appear to encounter values of 
7 where Tl is nondifferentiable in practice. 

We outline a scheme for calculating gradients of Ty and 
Tl in Appendix [A] The technique is essentially to use the 
implicit function theorem to find sensitivities of the variables 
in the problems that define Ty and Tl to the parameters 
7 , around the current solution of these problems, and then 
proceed to find the sensitivities of the optimal objective value 
for these problems to 7 . 
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(a) Even distribution of relative impedance changes 
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30 Bus / 41 Lines / 7 = 3.00 



10 20 30 40 

Number of Lines Attacked: k, 


(b) “Safe” distribution of relative impedance changes, to minimize total 
power adjustment 

Fig. 1. Power adjustment as a result of changes to line impedances on the 
30-Bus data set. It is often the case that load shedding is not required even 
when the disturbance to the system is quite substantial. 


C. Power-Adjustment Model Initialization 

As mentioned above, the objective value of the bilevel 
formulation is zero for most feasible values of 7 . It tends 
to be nonzero only on parts of the feasible region defined by 
(|1 lb[ ), ([TTc]) that correspond to near-maximal attacks focused 
on small numbers of buses. 

To illustrate this point, we perform experiments on the 30- 
Bus case (case30 . m from MatPower, originally from ESI) 
in which we monitor the power-adjustment objective Tl(j) 
in <[7]» as the impedances are increased. In Figure [T(a)| we plot 
Tl{i) in for the values 7 = 7 e, where 7 is a nonnegative 
scalar parameter that is increased progressively from 0 to 1. 
That is, impedances are increase evenly across all transmission 
lines. Note that T k ('ye) is zero for 7 £ [0,0.65], while for 
7 > 0.65, load shedding occurs on one or two demand buses. 
This observation implies that any value 7 along the line 7 e (for 
7 £ [0, .65]) is a global minimizer of the bilevel problem (|9](. 
The gradient is zero at each of these points, so optimization 
methods that construct the search direction from gradients 
cannot make progress if started anywhere along this line (or 
indeed from anywhere in a large neighborhood of this line). 
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If we are allowed to distribute a “budget” of impedance 
increases unequally between lines, so as to minimize the total 
amount of power adjustment required, even greater distur¬ 
bances can be tolerated. To describe this greater tolerance, we 
consider the following problem that is sliglty modified from 
©: 


Lower Bound Relaxed at Node 8 


■Hs(k, 7) : = min 

x,y, 7 

subject to 


T 

P V 

Fl(x,v; 7 ) = 0 
x<x<x 
0 < y < y 

T — 

e 7 = K 7 
0 < 7 < 7 e. 


(15a) 

(15b) 

(15c) 

(15d) 

05e) 

(15f) 


Note that © is different from 0 in two respects: (a) it 
is a single-level minimization problem whose variables are 
x , y. and 7; and (b) the budget is enforced with the equality 
constraint ( fl5el >. Thus this problem finds a “safe” way to 
distribute the fixed budget (^7) to transmission lines while 
the total load-shedding p T y is minimized. We solved this 
problem for an upper bound 7 = 3 with k, which is increased 
progressively from 1 to 41. The top chart in Figure |TTbJ| shows 
that it is possible to increase k to about 33 before any load 
shedding takes place at all. The lower chart depicts how the 
impedance changes are distributed around the 41 lines in the 
grid, at the solution of ( p~5j ), for each value of k. Darker bars 
on the graph show lines that can tolerate only a relatively 
small increase in impedance before causing load shedding 
somewhere in the grid. The lighter bars are those that can 
tolerate their impedance value 7 i being set to a value at or 
near the upper bound 3 without affecting load shedding. As 
an example: When n = 40, we have T~Ls( 40,3) ss 55, and the 
7 that achieves this power-adjustment value has components 
of 3 on all lines except line 16, where it is zero. 

The methodology used to derive Figure |l(b) can be used 
as a heuristic to identify a set of “safe” lines S (whose 
impedances can be increased without affecting grid perfor¬ 
mance) and a complementary set of “vulnerable” lines VV 
(for which impedance increases are likely to lead to load 
shedding). In Appendix [B] we describe the ESL (“estimating 
safe lines”) procedure, Algorithm[3j for determining the sets S 
and W. Once we have determined the vulnerable lines W, we 
define an impedance perturbation vector 7' with the following 
components: 

, i£VV 

7 , = < , (16) 

1 \o ii w, 

where 7 is the given upper bound on impedance on a given 
line. We then evaluate Fl{')') from 0- If a node does not 
require any load shedding under this maximal-perturbation 
setting, it is unlikely that any attack on the vulnerable lines 
will lead to load shedding on this node. We gather the other 
nodes — those for which /), > 0 at the solution of ( |7h| with 
7 = 7' — into a set T, the “target nodes.” Further explanation 
of the definition of T is given in Appendix |B2| 

We use the target nodes to define a modification of the 
objective F r that has the effect of shifting the range of the 



Lower Bound Relaxed at Node 20 


-10 


-20 



7io 


0 0 

(b) Non-target node 


740 


Fig. 2. Extending the range of 7) by allowing additional load on 
demand nodes. Adding load to target nodes (top figure) produces a useful 
extension of the range of the objective, so that its gradient yields a promising 
search direction for the maximization problem. Adding load to non-target 
nodes (bottom figure) simply shifts the objective down by a constant, so that 
the gradient in the flat region still does not yield a useful search direction for 
the maximization problem. 


function, in a way that makes gradient information relevant 
even at values of 7 for which no power adjustments are 
required. The idea is illustrated in Figure [2j where we show the 
power-adjustment requirement on two different nodes (nodes 
8 and 20) of the 30-Bus system as a function of the values of 
two impedance parameters — those corresponding to lines 10 
and 40. In both graphs of Figure [2] the top surfaces (shaded 
white and red) represent the objective Fl as a function of 
various values of the pair (710,740)- Note that Fl takes the 
value zero over much of the domain, but becomes positive 
when both 710 and 740 are high. The lower surfaces in each 
graph show how Fl changes when we modify the subproblem 
in 0 by removing the zero lower lower bound on the load 
Pi in 0t]i, where i = 8 (a target node) in Figure |2(a)| and 
i = 20 (a non-target node) in Figure |2(b) Removal of the 
lower bound has the effect of allowing load to be added to 
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Algorithm 2 Vulnerability Analysis: Power- 
Adjustment Model_ 

Require: 

7 : Upper bound for impedance increases 7 $, i E £; 

Number of lines to attack; 

7 0 : Feasible initial value of 7 ; 

Ensure: 

7 *: Impedance vector that optimizes the attack; 


iPOPT^CWachter and Biegler H 1 7| ) as the nonlinear solver for 
evaluating Fl 0 in the power-adjustment (bilevel) model. 
For the test case data and calculation of the electric circuit 
parameters, the codes from MatPower lfl8l are used exten¬ 
sively. The codes were executed on a Macbook Pro (2 GHz 
Intel Core i7 processor) with 8GB RAM. 


1: Find a set of vulnerable lines W and target nodes T using the ESL 
procedure (Algorithm [3j; 

2: Set lower bound of pi (for target nodes i E 1~) to — oo; 

3: k «- 0; 

4: while k <MaxIter do 
5: Find gradient g k of Fl at 7 . 

6 : Find linearized optimum w k from { 0 ; 

7: Use the backtracking to find step size a^ E [0,1]; 

8: y=+i ^ 7 fc + Qfc( y=_ 7 7 

9: Identify i such that pi = argminpy, where pj are the 

i 

power-adjustment variables from (7); 

10: if pi < 0 then 

11 : reset lower bound on pi to zero; 

12: end if 

13: ki-k + 1: 

14: Stop if terminating conditions (including pj > 0 for all power- 

adjustment variables pj) are satisfied; 

15: end while 

16: 7 * «— 7 fc ; 


the node in question. This is not an action that would be 
operationally desirable, but as we see from the blue surface 
in Figure 2(a) it changes the nature of Fl in useful ways. 
The effect of removing the lower bound on p 8 (Figure 2(a) i 
is to extend the range of Fl so that its derivative at any point 
in the domain gives useful information about a good search 
direction. In a sense, the extended-range version appears to be 
a natural extension of the original objective Fl- By contrast, 
removal of the lower bound on p 2 o (Figure |2(b)[ > causes the 
function to simply be shifted downward by a roughly constant 
amount for all pairs of impedance perturbation values. This is 
because, being a non-target node, increased load on this node 
can be met, even after the grid is damaged by the impedance 
attack. We conclude that removing lower bounds on pi for 
target nodes i G T provides a potentially useful extension of 
the range of the function Fl, whereas the same cannot be said 
for non-target nodes. 

Motivated by these observations, we modify Algorithm |T] 
as follows. We start by removing all lower bounds in ( |7h| on 
target nodes i G T- At each iteration of the algorithm, after 
taking a step, we check to see if any of the pi obtained by 
solving the subproblem (|7]» at the latest iteration are negative. 
If so, we reset the lower bound on the most negative value of pi 
to zero, before moving on to the next iteration. The algorithm 
does not terminate until all p, are nonnegative. The modified 
procedure is shown as Algorithm [2] 


IV. Experimental Results 

We present the results obtained with our formulations 
and algorithms on the IEEE 118-Bus system and Polish 
2383-Bus system. Our implementations use MatlaeQ with 


TABLE I 

Test Cases for Experiments 



Test Cases 

1 1 2 

3 | 4 

Filename (in MatPower) 

case 118.m 

case2383wp.m 

Number of Nodes 

118 

2383 

Number of Lines 

186 

2896 

Number of Lines to Attack k, 

3 | 5 

3 | 5 

Perturbation Upper Bound 7 

3 

2 

Backtracking Parameters (co, ci, o; m i n ) 

(0.5, 0.01, 0 . 01 ) 

Voltage limits (V, V) 

(0.93,1.07) | (0.89,1.12) 

Line Screening Threshold 77 

0.9 


Information about the test case instances and algorithmic 
parameters are given in Table [I] There are two instances for 
each of the two grids, corresponding to 3-line and 5-line 
attacks, respectively. The table shows voltage magnitude limits 
that are applied in the power-adjustment model, together with 
the value of rj that is used in the ESL procedure (Algorithm [3] 
from Appendix |B1| >. 

In the power-adjustment model |7]), the upper bounds of, 
of, and Pi on the power-adjustment variables are set to 1 for 
most buses, thus allowing full load shedding. If a bus violates 
our assumption on power injection — that is, if Pi < 0 for 
bus i G G or Pi > 0 or Qi > 0 for bus i G V — the load¬ 
shedding upper bound for that bus is set to 0, disallowing 
power adjustment on that bus. The power-adjustment objective 
Fl is considered to be nonzero if it is at least 10” 3 megawatt 
(MW). 


A. Voltage Disturbance Model 

We discuss first results obtained with the voltage disturbance 
model <[4}-([5]i applied to the four test cases of Table [I] 

118-Bus System: For a 3-line attack problem on IEEE 118- 
Bus system (k = 3), Algorithm [I] converges in 5 iterations and 
identifies exactly three lines to attack with maximal impedance 
increase: lines 71, 74, and 82 (as shown in Table [II]i. As a 
result of this attack, voltage magnitudes at four buses decrease 
significantly, by up to 0.07 p.u., as shown in Figure [3] (In 
Figure |3(a) the buses are reordered in increasing order of 
voltage magnitude on the undisturbed system. In Figure |3(b)| 
the buses are indexed in their original order.) The attack is 
visualized in Figure[5] where we see that its effect is essentially 
to isolate buses 51, 52, 53, and 58; the attacked lines are 
colored in red. 

For the second test instance, on the IEEE 118-Bus system 
with k = 5, the algorithm identifies exactly five lines to attack 
at the maximal impedance increase — lines 25, 29, 71, 74, and 
82 (see Table — which includes the three lines identified 
in the first attack. With this stronger attack, there is significant 
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TABLE II 

Voltage Disturbance Model: 118-Bus System with k = 3 


TABLE III 

Voltage Disturbance Model: 118-Bus System with k = 5 


Line 

No. 

Buses 

Continuous 
Attack (7 i) 

From 

To 

71 

49 

51 

3.00 

74 

53 

54 

3.00 

82 

56 

58 

3.00 

Objective 

4.04 x 10~^ 


Optimal Attack (as determined by our algorithm) 
Voltage Magnitude 



x After Attack 
• Before Attack 


20 40 60 80 100 

Bus Index (reordered) 

(a) Distribution of Voltage Magnitudes Before and After Attack. 

Changes in Voltage Magnitude 


Line 

No. 

Buses 

Continuous 
Attack (7 i) 

From 

To 

25 

19 

20 

3.00 

29 

22 

23 

3.00 

71 

49 

51 

3.00 

74 

53 

54 

3.00 

82 

56 

58 

3.00 

Objective 

5.03 x 10~' 2 


Optimal Attack (as determined by our Algorithm) 
Voltage Magnitude 



x After Attack 
• Before Attack 


20 40 60 80 100 

Bus Index (reordered) 

(a) Distribution of Voltage Magnitude Before and After Attack. 


0.00 - - 

- 0.02 

3 

Cl, 

-0.04 

-0.06 * * 

_I_I_ •! _I_I_1_ 

0 20 40 60 80 100 

Bus Index (original order) 


(b) Distribution of Voltage Magnitudes Changes After Attack. 
Fig. 3. Voltage Disturbance Model: 118-Bus System with k, = 3 


voltage drop on seven buses. Algorithm [I] takes 8 iterations 
to converge to the solution. As Figure [5] shows, attacking 
the additional two lines (colored in green) has the effect of 
creating another “island,” consisting of buses 20 , 21 and 22 . 
The additional voltage drops seen Figure |4(b) are from these 
buses. 

2383-Bus System: Results for our third test instance in 
Table H] which considers the 2383-Bus model with attack limit 
defined by k = 3, are shown in Table IV The continuous 
impedance attack is distributed into 5 lines, identified after 
ten iterations of Algorithm |T| (The lines involved in the attack 
are revealed at iteration five, while the remaining five iterations 
make minor adjustments to the impedance values.) 

This 5-line solution can be used to identify the most 
disruptive set of three lines using two heuristics: (a) choose 
the three lines i such that 7 ,; are one of the largest three entries 
in 7 (called the Top-3 attack ); and (b) try all possible 3-line 
combinations of the five lines with highest impedance (the 
Best-3 attack). For this specific case, the Top-3 attack and 
Best-3 attack are the same, consisting of lines 5, 405, and 467. 
These 3-line attacks, with impedance set to their maximum 
values on all three lines, gives a slightly smaller objective 
value than the continuous attack. 


Changes in Voltage Magnitude 


0.00 

- 0.02 
3 

Cl, 

-0.04 
-0.06 

0 20 40 60 80 100 

Bus Index (original order) 



(b) Distribution of Voltage Magnitude Changes and After Attack. 
Fig. 4. Voltage Disturbance Model: 118-Bus System with k = 5 


Figure | 6 (bj| shows how the voltage magnitude changes when 
the Best-3 and Top-3 attacks are executed. Similarly to 118- 
Bus cases, there is a relatively small number of buses in which 
the voltage drops significantly from its original value, but large 
voltage changes are seen on some buses, with some voltage 
magnitudes below 0.8 p.u. 

For our fourth test instance in Table [I] for which k = 5 
on the 2383-Bus system, two iterations of Algorithm [l] suffice 
to identify an attack (on lines 404, 405, 467, 479, and 501) 
that makes the power flow problem infeasible, that is, there 
is no (V,6) that satisfies F(V,9; 7 ) = 0 under this attack. 
Hence, unless the grid operator takes action (to change loads 
or generator outputs, for example), an attack on these five lines 
renders the grid inoperable. 


Comparison with Continuous Optimization Model: In Sec- 
tion |II-B we mentioned that the voltage disturbance model can 
be written as a continuous optimization model ([ 6 |, except that 
the latter model does not handle infeasibility appropriately. We 
verify the properties of this alternative formulation by solving 
it with the nonlinear interior-point solver Ipopt. We find that 
in the first three instances of Table |IJ the solutions obtained 
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Fig. 5. Voltage Disturbance Model: Attacks on the 118-Bus System. The transmission lines in red are are the optimal attack for k = 3. The lines in green 
are added to the optimal 3-line attack when k is increased to 5. 


TABLE IV 

Voltage Disturbance Model: 2383-Bus System with k = 3 


Line 

No. 

Buses 

Continuous 
Attack ( 7 i) 

Top-3 

Attack 

Best-3 

Attack 

From 

To 

5 

10 

3 

1.04 

2.00 

2.00 

404 

434 

188 

0.25 



405 

437 

188 

2.00 

2.00 

2.00 

467 

340 

218 

2.00 

2.00 

2.00 

501 

340 

240 

0.71 



Objective 

0.514 

0.501 

0.501 


Optimal Attack (as determined by our Algorithm) 


Voltage Magnitude 


1.1 



0.9 L “ 


J 


x After Attack 
• Before Attack 

_J_l_I_ . 

500 1,000 1,500 2,000 

Bus Index (reordered) 

(a) Distribution of Voltage Magnitude Before and After Attack. (Best-3) 



Changes in Voltage Magnitude 

0.00 ■ 

. -0.05- 7 

3 

Cl. 

- 0.10 

-0.15 |~ '* * ■>*. _|_i_ i 

0 500 1,000 1,500 2,000 

Bus Index (original order) 

(b) Distribution of Voltage Magnitude Changes After Attack. (Best-3) 
Fig. 6 . Voltage Disturbance Model: 2383-Bus System with k, = 3 


from © match those we described above. In the fourth test 
case — the 2383-Bus model with k = 5 — the model ([6]) 
identifies a solution that maximizes disruption subject to the 
power flow equations {3} remaining feasible. Our model ([5}, 
which detects infeasibility of the grid under an attack of this 
strength, yields the more informative outcome. 


B. Power-Adjustment Model 

We now present results for the power-adjustment model |7|, 
for the four test instances of Table [I] 

118-Bus System: The ESL procedure (Algorithm [3]), which 
is described in Section III-C and Appendix |B1| is applied 
to the 118-Bus system to identify vulnerable lines VV = 
{29, 71, 74,96,184} and target nodes T = {52, 53,117}. We 
use these settings of VV and T in Algorithm [2] to solve the 
first two instances in Table [j] 

For the first test case (k = 3), Algorithm [2] terminates after 
eight iterations, at the attack shown in Table [V] Four lines 
are involved in this attack; we reduce to three-line attacks 
“Top-3” and “Best-3” as in Subsection IV-A The Top-3 and 
Best-3 attacks coincide (and are the same as those obtained for 
the voltage disturbance model in Subsection IV-Ai and have 
a slightly smaller objective value than the continuous attack. 
We note that in these 3-line attacks, the voltage magnitudes of 
buses 52 and 53 are at their lower bound V = .93, and load 
shedding is required for buses 51 and 53, the total amount of 
load shedding being 22.13 MW. 

In Table VI we show the top five three-line attacks obtained 
by enumerating all 266,916 possible three-line attacks on this 
grid. The most disruptive attack is indeed the one found by 
our algorithm. 

For the second test instance (with k = 5), Algorithm |T] 
requires twelve iterations to converge, and distributes the 
impedance increases around nine lines, with the maximum 
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TABLE V 

Power-Adjustment Model: 1 18-Bus System with k = 3 


TABLE VIII 

Power-Adjustment Model: 2383-Bus System with k = 3 


Line 

Buses 

Continuous 

Top-3 

Best-3 

No. 

From 

To 

Attack ( 7 i) 

Attack 

Attack 

71 

49 

51 

3.00 

3.00 

3.00 

72 

51 

52 

0.26 



74 

53 

54 

3.00 

3.00 

3.00 

82 

56 

58 

2.74 

3.00 

3.00 

Objective (MW) 

22.39 

22.13 

22.13 

Buses with 

Load Shedding 

51, 52, 53 

51, 53 

51, 53 

Buses at 

Voltage Boundary 

52 

52, 53 

52, 53 


TABLE VI 

Five Most Effective Attacks from Exhaustive Enumeration: 
1 18-Bus System with k = 3 


Lines Selected 

Power Adjustment (MW) 

71 

74 

82 

22.13 

71 

72 

74 

19.07 

71 

74 

83 

17.21 

71 

74 

184 

15.61 

71 

74 

97 

13.27 


perturbation (7,; = 3) on four of them; see Table VII The Top- 
5 and Best-5 attacks each involve the four lines with maximum 
impedance increases, but differ in their choice of additional 
line. Three buses require load shedding in the Best-5 attack, 
compared to two in the Best-3 attack. 

Unlike the previous results for k = 3, which target the 
same lines as the voltage disturbance model, the Best-5 attack 
identified in the power-adjustment model is slightly different 
from the Best-5 attack for the voltage-disturbance model. The 
power-adjustment attack targets the region around buses 51, 
52, and 53, as before, but also “islands” bus 117 (see Figure[5]l, 
rather than the buses 20, 21, and 22 that are attacked by the 
voltage-disturbance model. 

2383-Bus System: For the Polish 2383-Bus system, the ESL 
procedure (Algorithm [3} identified a set W of 20 vulnerable 
lines and, for upper bound 7 = 2 on the impedance increase, 
a set T of 55 target nodes. (Algorithm [3] required about 270 


TABLE VII 

Power-Adjustment Model: 118 -Bus System with k = 5 


Line 

Buses 

Continuous 

Top-5 

Best-5 

No. 

From 

To 

Attack ( 7 i) 

Attack 

Attack 

71 

49 

51 

3.00 

3.00 

3.00 

72 

51 

52 

0.36 


3.00 

74 

53 

54 

3.00 

3.00 

3.00 

75 

49 

54 

0.74 



76 

49 

54 

0.81 

3.00 


81 

50 

57 

0.29 



82 

56 

58 

3.00 

3.00 

3.00 

97 

64 

65 

0.80 



184 

12 

117 

3.00 

3.00 

3.00 

Objective (MW) 

27.16 

25.79 

26.87 

Buses with 
Power Adjustment 

51, 52, 53, 117 

51, 53, 117 

52, 53, 117 

Buses at 

Voltage Boundary 

52, 117 

52, 53, 117 

52, 117 


Line 

Buses 

From Bilevel Formulation 

From N - 1 

No. 

From 

To 

Cont. ( 7 i) 

Top-3 

Best-3 

Top-3 

Best-3 

5 

10 

3 




2.00 


264 

140 

117 

0.02 





268 

126 

118 

2.00 

2.00 

2.00 

2.00 

2.00 

289 

135 

125 

1.98 

2.00 

2.00 


2.00 

296 

145 

128 

2.00 

2.00 

2.00 

2.00 

2.00 

Objective (MW) 

577.80 

577.75 

577.75 

303.07 

577.75 

# of Buses with 
Power Adjustment 

49 

27 

49 

Buses at 

145, 

146, 1905 

145, 146, 

145, 146, 

Voltage Boundary 

230, 1905 

1905 


seconds of run time on this instance.) 

Results for the power-adjustment model on or our third test 
case from Table [Tj for the 2383-Bus system with k = 3, are 
shown in Table VIII Algorithm [2] returns a four-line attack. 
Since the attack on one of these four lines (line 264) is 
negligible, we find that the Top-3 and Best-3 solutions both 
attack lines 268, 289 and 296, with an active-power load 
shedding 577.75 MW, which is negligibly smaller than the 
optimal four-line attack. There are 49 buses which need load 
shedding for the three-line attack, and three buses have voltage 
magnitudes at their lower limits of V = .89 — another sign 
of stress on the grid. The solution obtained for this test case 
is quite different from the one from the voltage disturbance 
model. The lines 5, 404, 405, 467, and 501 which are identified 
by the voltage disturbance model (cf. Table IV 1 also cause 
some load shedding when they are attacked, but the effect is 
not as serious as the attack on lines 268, 289, and 296. 

Since it is computationally intractable to look at all possible 
three-line combinations in a 2383-Bus grid, we apply an 
“N — 1 enumeration” heuristic to explore the most promising 
part of the space of three-line attacks. In this heuristic, each 
line i is individually perturbed by setting 7,; = 7, and we 
note which of these perturbations require load shedding. On 
this data set, sixteen lines were identified as causing load 
shedding. We define a “Top-3 (N — 1)” attack to comprise 
the three lines that individually cause the most load shedding, 
and the “Best-3 (AT — 1)” attack to be the most dismptive 
three-line combination drawn from these sixteen lines. The 
resulting attacks are displayed in Table VIII alongside the 
Top-3 attack and Best-3 attack. We see that the Best-3 (IV — 1) 
attack is identical to the Top-3 and Best-3 attacks, while the 
Top-3 (N — 1) attack is inferior. 

The optimal attacks for the power-adjustment model in this 
third test instance are quite different from those obtained 
from the voltage disturbance model, as we see by comparing 
Tables [IV] and |VI1I| We found, however, that if the lower 
bound on voltage magnitude is changed from V = .89 to 
V = .87, the optimal attack for the power-adjustment model 
is almost identical to the voltage disturbance model. In the 
relaxed problem, the buses 145, 146, and 1905 no longer have 
their voltage magnitude at the lower bound, while buses 401 
and 414 move to the relaxed lower bound. The latter two buses 
are among those that suffer significant voltage drop in the 
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TABLE IX 

Power-Adjustment Model: 2383-Bus System with k = 5 


Appendix 


Line 

Buses 

From Bilevel Formulation 

From N — 1 

No. 

From 

To 

Cont. (74) 

Top-5 

Best-5 

Top-5 

Best-5 

5 

10 

3 




2.00 


268 

126 

118 

2.00 

2.00 

2.00 

2.00 

2.00 

269 

142 

118 

2.00 

2.00 

2.00 



289 

135 

125 

2.00 

2.00 

2.00 

2.00 

2.00 

296 

145 

128 

2.00 

2.00 


2.00 

2.00 

317 

142 

135 

1.54 

2.00 

2.00 



405 

437 

188 




2.00 

2.00 

467 

340 

218 





2.00 

2142 

1693 

1658 

0.46 


2.00 



Objective (MW) 

1109.72 

1086.67 

1460.25 

594.09 

597.71 

# of Buses with 
Power Adjustment 

77 

78 

71 

53 

53 






145, 146, 

145, 146, 

Voltage Boundary 

145, 146, 1905 

1905 

230, 414, 
434, 1905 

401, 414, 
1905 


optimal voltage-disturbance attack of Table IV 


Table IIXl shows results for our fourth test instance from 
Table |T| an optimal attack on the 2383-Bus system for k = 5. 
The attack determined by our procedure is distributed around 
6 lines. The Top-5, Best-5, Top-5(./V — 1), and Best-5(7V— 1) 
attacks are calculated as described above. The Best-5 attack 
is significantly more disruptive than the optimal continuous 
attack identified by Algorithm [2] probably because of noncon¬ 
cavity in the objective Tl- We note however that our algorithm 
is much more useful in screening for the most disruptive 
collection of lines than is the standard “N — 1” screening 
methodology: The Top-5 and Best-5 attacks are much more 
damaging than the Top-5 (N — 1) and Best-5 (N — 1) attacks. 


V. Conclusions 

We have proposed an attack model for assessing the vul¬ 
nerability of power grids. The attack consists in increasing the 
impedance of transmission lines, with the resulting disruption 
to the grid measured in two ways. The first technique is to 
observe changes in voltage magnitudes at the buses; greater 
changes from the nominal values indicate greater disruption. 
The second technique is to measure the weighted sum of 
adjustments to load and generation that are needed to restore 
stable operation of the grid, with the voltage magnitudes 
confined to certain prespecified ranges. The two criteria give 
rise to optimization problems with different properties. Both 
are solved with a combination of known algorithms (such as 
Frank-Wolfe) and heuristics that determine promising regions 
of the solution space. 

In our computational results, we also use our algorithm as a 
screening procedure for determining which collections of lines 
are likely to cause the most disruptive attacks. By enumerating 
combinations of lines from among those identified by our 
algorithms, we identify more disruptive attacks than those 
produced by alternative screening methods, such as the well 
known “N — 1” criterion. 


A. Computing Gradients 

We describe here calculation of gradients for the functions 
T ( 7 ) defined in Section [IT] which quantify the grid disruption 
arising from an attack modeled by the impedance vector 7 . 
Both model functions considered here — 0 and (|7]i — have 
the following form: 

F{l) := min f(z) (17a) 

Z 

subject to Ci(z, 7 ) = 0, i = l,2,...,m, (17b) 

hj(z)> 0, j = 1,2,... ,r, (17c) 

where the functions /, < 7 , i = 1 , 2 and hj, j = 
1,2 ,...,r are all smooth. (Note that in our models, the 
dependence on the upper-level variables 7 arises only through 
the equality constraints, but our discussion can be extended 
without conceptual difficulty to the case in which the objective 
and the inequality constraints also depend on 7 .) 

We outline a technique for calculating the gradient VJ r ( 7 ). 
We assume that the minimizing z for the (generally noncon- 
vex) problem has been identified and that it is denoted 
by z( 7 ). Moreover, we assume that z( 7 ) is a nondegenerate 
solution of the minimization problem in G3- By this we 
mean that the linear independence constraint qualification 
holds at the minimizer, that a strict complementarity condition 
holds, and that second-order sufficient conditions are satisfied 
at z( 7 ). We show that under these conditions, we can use 
the implicit function theorem to define the gradient VJ r ( 7 ) 
uniquely. While strong, these conditions are not impractical; 
they appear to hold for all values of z( 7 ) encountered by our 
algorithm, and it seems plausible that they would hold for 
“almost all” values of 7 . The question of existence of VJ r ( 7 ) 
becomes much more complicated when these conditions are 
not satisfied. When strict complementarity does not hold, for 
example, the set of active inequality constraints is on the 
verge of changing, an event usually associated with a point 
of nonsmoothness of 

The Karush-Kuhn-Tucker (KKT) conditions for optimality 
of z in the problem ( fTT) for a given 7 are that there exist 
scalars A*, i = 1 , 2 ,..., to and pj, j = 1 , 2 ,..., r such that 

m r 

V/(z) - X>V*Ci(*, 7 ) - = 0. ( 18a ) 

*=1 i=1 

Ci(z, 7 ) = 0, i = 1,2,... ,m, (18b) 

Pj > 0 , j e A(z), (18c) 

hj{z) = 0 , jeA{z), (18d) 

Pj = 0, j G {1,2, ...,r}\A(z), (18e) 

hj(z)> 0, j G {1,2,...,r} \ A{z), (18f) 

where the active set A(z) is defined as follows: 

A(z) :={j = l,2,...,r : h 3 {z) = 0}. (19) 

We use the following vector notation: 

c(z, 7 ) = [a(z, 7 )]™!, A = [A,]™ 1; 

dA = [dj}jGA, h A {z) = [hj(z)\ jeA . 
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The linear independence constraint qualification (LICQ) is: 

{V~Ci(z, 7 ), i = 1,2,.. .,m}U{V/ij( 2 ), j G -4(z)} (2() , 
is a linearly independent set. 

The strict complementarity condition is that 

p : j > 0, for all j G A{z). (21) 

Finally, the second-order sufficient conditions are 

d T W(z, X, p, 7 )d > 0 for all d G N(z) with d ^ 0, (22) 

where the subspace N(z) is defined as follows: 

f V z Ci{z, 7 ) T d = 0 , i = 1 , 2 ,..., TO 1 

iV(«) := < d : and J , (23) 

( Vhj{z) T d = 0, Vj G .A(«) J 

and the matrix VF (z, A, p, 7 ) is the Hessian of the Lagrangian 
function for the problem (jT7j, that is, 


W(z,X,p, 7) := 

m 

i=i ie.4(z) 


(24) 


In the neighborhood of a value of 7 at which all the 
conditions above are satisfied, we find that (z, A, p A ) is an 
implicit function of 7. We find expressions for the derivatives 
of (z,X,p A ) with respect to 7 by applying the implicit 
function theorem (see for example 119] Theorem A.2]) to the 
equality conditions in (fl8|), which we can formulate as follows: 


v /(-)^E AiVzCi ^ ,7 ) _ E t jL j^ h j( z ) = 0 ’ (25a ) 

*= 1 jeA(z) 

Ci(z,nr) = 0, i = 1,2,... ,m, (25b) 

hj(z)= 0, j G A(z). (25c) 

Note that this system of equations is square, with n + m + 
\A\ equations and unknowns. Moreover, standard analysis of 
optimality conditions shows that its square Jacobian matrix is 
nonsingular, under the LICQ ( |20| i and second-order sufficient 
( |22| > conditions. The implicit function theorem now yields the 
following: 


(26) 


’ Vz( 7) ' 


~ m 

^ ^ Aj V z'fC'i (z 5 t) 

VA( 7 ) 

= H(z, X, p, 7) 1 

i =1 

T—7 / \T 

y pa ( 7 ). 


—V 7 c(z, 7 ) t 

0 


where 


H(z,\,p,j) = 


W{z,X,p, 7 ) -V z c(z, 7) -Vh A {z) 
X7 Z c(z, 7 ) t 0 0 

S7h A {z) T 0 0 


We can derive VJ r (z) from Vz( 7 ) through the definition ( fl7| , 
as follows: 

VT( 7 ) = X7z( 7 ) T Vf(z( 7 )). (27) 


B. Heuristics for Initializing the Power Adjustment Model 

1) Determining Safe and Vulnerable Lines: We have noted 
that a typical grid contains many “safe” lines, for which large 
changes to the impedance do not affect the ability of the grid 
to serve demands. We discuss here a filtering approach to 
identify the complementary set of “vulnerable” lines, for which 
impedance change causes significant disruption of the grid. 
(We assume that the number of vulnerable lines is relatively 
small; otherwise, the grid has a systemic vulnerability and 
would be hard to defend.) 

A naive approach for identifying safe and vulnerable lines 
is enumeration. For each line i, we set 7 = 7 e t (where e, is a 
vector of all zeros except for 1 in the ith entry) and evaluate 
(Fl{ 7 ) defined by ([8]). Vulnerable lines are taken to be those 
for which Fl{ 7 ) > 0. This enumeration approach is not very 
effective, in part because of its cost (it requires solution of £ 
different power flow problems) and because it cannot identify 
combinations of lines that are individually “safe” but which 
together create a vulnerability in the network. We therefore 
propose an alternative heuristic called ESL (for “eliminating 
safe lines”). 

The motivation of the ESL heuristic is as follow. If a system 
operator, instead of an attacker, is asked to increase impedance 
on exactly k lines, then he will choose those k lines so as 
to minimize the disruption to the system. These lines will 
not be an attractive choice on the attacker’s side, especially 
when increase of impedances on these lines does not lead to 
any load shedding, so in general we declare these lines to 
be “safe.” On the other hand, the lines not selected by the 
operator correspond to those that may cause disruption and 
hence are an attractive target for attack. We declare such lines 
to be “vulnerable.” 


In ESL, following the experiment graphed in Figure 1(b) 


we seek the value of n in ( |9b| ) such that a total perturbation 
of size n 7 can be distributed to lines with little load shedding. 
Lines i for which 7 , ft: 7 are declared to be safe. This process 
is repeated until relatively few vulnerable lines remain. The 
appropriate value of k can be found by binary search, by 
solving the following modification of the problem in 
which depends on a working set W C C of lines not yet 
classified as safe: 


Kw(K)7) = nrin 

x,y, 7 

T 

v y 


(28a) 

subject to 

11 

'p 

a 

0 

(28b) 


e T 7 = n 7 


(28c) 


x<x<x 


(28d) 


0 <y<y 


(28e) 


0 < 7» < 7 

i G W 

(28f) 


7 i = 0 

i i w. 

(28g) 

The complete procedure is shown in Algorithm [3] We start 
by putting all lines £ into the working set W, then successively 


eliminating from W those lines i for which the solution of ( 
yields 7 * within a factor // of the upper bound 7 . (We used 
?/ = .9.) The process is repeated until no new “safe” lines 
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Algorithm 3 Eliminating Safe Lines (ESL) 

Require: 

£ : Set of all lines; 

r\ S (0,1): Threshold for screening; 

Ensure: 

W: a set of vulnerable lines; 

S'. a set of safe lines; 

1 : 5 <— 0 ; 

2 : repeat 

3: W i- C\S; 

4: Define k* to be the largest value of k for which Pyv( K > 7) < e; 

5: S' 3— {i : 7^/7 > r], i G W}; ► newly determined “safe” lines 

6: S^SUS'; 

7: until S 1 = 0 



Fig. 7. Loadability and Load Shedding in PV-curve 



(a) Non-Target Node (No Load Shedding Required) 



(b) Target Node (Load Shedding Required) 


Fig. 8. Possible Changes of PV-curve on a Node After Attack. 


are identified. The lines remaining in W are then classified as 
“vulnerable.” 

2) Target Node Selection: Our approach for selecting “tar¬ 
get” nodes T in Algorithm [ 2 ] is based on maximum loadability. 
The maximum loadability problem is similar to feasibility 
restoration in that it seeks the boundary of the feasible region. 
However, rather than starting from an infeasible point (where 
the nominal loads cannot be served), it begins from a feasible 
grid and increases the loads until demands can no longer be 
met. The difference is illustrated in Figure [7] which shows 
the PV-curve for a particular demand node. When the grid 
is feasible with demand Pd, (right curve), the demand can 
be increased to P'/j while retaining feasibility. The difference 
P'd — Pd can be regarded as the maximum loadability at this 
node. If the grid is infeasible (left curve) the demand must be 
reduced to P' D before feasibility is recovered. 

The formulation for maximum loadability problem can be 
obtained by replacing ( |7f| ), ( |7g[ >. and ( |7h| by the following 
constraints: 


a+= 0 

i G Q 

(29a) 

G~ = 0 

i G Q 

(29b) 

o 

VI 

i G V. 

(29c) 


The first two constraints fix the power generations at their 
nominal values, while \29c) allows increase (rather than de¬ 
crease) of demand at the demand buses. When the nominal 
loads and generations are feasible, we expect the objective to 
be negative at the solution. 

To identify the target nodes, we simply set 7 i = 7 for 
all vulnerable lines i G W that are identified by the ESL 
procedure. Algorithm [3] and solve 0 for this value of 7. 
If a node does not require any load shedding under this 


maximal-perturbation setting, it is unlikely that any attack 
on the vulnerable lines will lead to load shedding on this 
node. The target nodes are defined to be those for which load 
shedding is required, that is, pi > 0 at the solution of ([TJ). We 
denote the set of these nodes by T. 

Figure [ 8 ] shows target and non-target nodes, and shows 
how maximum loadability motivates their classification into 
these categories. The top figure shows a non-target node, for 
which is it possible to meet the original demand Pd even 
after the maximum-perturbation attack, though the loadability 
is decreased. The bottom figure shows that the demand must 
be reduced to P'd in order for the network to remain feasible. 
On this node, there is a chance that an attack on the vulnerable 
lines will lead to load shedding. The target nodes are the nodes 
affected by the type of attack we are considering, so that even 
for a choice of 7 that allows the nominal demand on these 
loads to be served, the change in maximum loadability may 
give us some information on the sensitivity of demand that 
can be served to the value of 7. Since the objective in ([7]) is 
of weighted £1 type, we expect the number of target nodes to 
be small. 

The target nodes T are incorporated into Algorithm [2] 
by replacing the lower bound © in the formulation 0 ) 
by a negative quantity for the nodes in T, and increasing 
these bounds toward zero progressively during the course of 
Algorithm [2] Additional details are given in Subsection III-C 
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